close all
clear all
Fig_settings
% GRID PARAMETERS
n=160;
N=n;
R=2.0;
mu=1;
epsilon=0.25*2*pi*R/N;
t_step=5;
% BOUNDARY PARAMETERS
t=[0:2*pi/(t_step-1):2*pi];
%t=0;
om=1;
% IMPOSED VELOCITY
u=1.0;
v=0.0;
for j=1:N
    vel_in((j-1)*2+1)=u;
	vel_in((j-1)*2+2)=v;
end
% GRID POINTS ON A CIRCLE
[xb,yb,theta,dtheta]=long_lat_coord2D(n,R);
 % BREATHING BOUNDARY
 
 for i=1:t_step
     
     [xb,yb,ub,vb,theta,dtheta]=moving_coord2D(n,R,t(i),om);
     %clf;
     figure
     plot(xb, yb, 'k.')
     title(['t= ' num2str(t(i))])
     xlabel('X')
     ylabel('Y')
     daspect([1 1 1])
     axis([-1.5*R 1.5*R -1.5*R 1.5*R])
     grid on
    % pause(.1)
    % A(:,i) = getframe(gcf);
    % movie2avi(A, 'boundary.avi', 'fps', 5, 'quality', 100);
% THE REGULARIZED GREEN FUNCTION MATRIX IS COMPUTED
[G]=G_matrix_mod2D(xb,yb,N,xb,yb,N,mu,epsilon);

figure
imagesc(G)
colorbar
title('G')

% SOLVE THE LINEAR SYSTEM TO FIND THE FORCING f ON THE BOUNDARY
opts.SYM = true;
f=linsolve(G,vel_in',opts);
 
% FIND VELOCITIES IN THE FLOW DOMAIN ONCE THE FORCING ON THE BOUNDARY IS
% KNOWN

%xout=R*[cos(pi/128),sin(pi/128); cos(pi/12),sin(pi/12)];
%xout=[2*R,0];
%[M2N]=G_matrix_2XN(xout,xb,yb,N,mu,epsilon);

%%% CALCULATE VELOCITY ON A GRID 
xout=[-3:0.25:3];
yout=[-3:0.25:3];
[xoutM, youtM]=meshgrid(xout,yout);
Nout=length(xout);
[M2N]=G_matrix_mod2D(xoutM(:),youtM(:),Nout^2,xb,yb,N,mu,epsilon);

%%% CALCULATE VELOCITY IN FEW POINTS
%xoutM=R/8*[cos(pi/3) cos(pi/6)];
%youtM=R/8*[sin(pi/3) sin(pi/6)];
%Nout=length(xoutM); 
%[M2N]=G_matrix_mod2D(xoutM,youtM,Nout,xb,yb,N,mu,epsilon);

uout=M2N*f;
%for j=1:N
%    f_an2((j-1)*2+1)=8*pi/(1-2*log(R))/N;
%	f_an2((j-1)*2+2)=0;
%end
%uout2=M2N*f_an2';
uout_vec=uout(1:2:2*Nout^2);
vout_vec=uout(2:2:2*Nout^2);
for i=1:Nout
    uoutM(1:Nout,i)=uout_vec((i-1)*Nout+1:i*Nout);
    voutM(1:Nout,i)=vout_vec((i-1)*Nout+1:i*Nout);
end

%uout_vec2=uout2(1:2:2*Nout^2);
%vout_vec2=uout2(2:2:2*Nout^2);
%for i=1:Nout
%    uoutM2(1:Nout,i)=uout_vec2((i-1)*Nout+1:i*Nout);
%    voutM2(1:Nout,i)=vout_vec2((i-1)*Nout+1:i*Nout);
%end
%%% ANALYTIC SOLUTION
f_an=8*pi/(1-2*log(R)).*[1,0];
modx2=xoutM.^2+youtM.^2;
uout_an=-f_an(1)/(8*pi)*(2.*log(sqrt(modx2))-R^2./(modx2))...
        +(f_an(1).*xoutM+f_an(2).*youtM).*xoutM./(4*pi.*modx2).*(1-R^2./modx2);
     
vout_an=-f_an(2)/(4*pi)*(2.*log(sqrt(modx2))-R^2./(modx2))...
        +(f_an(1).*xoutM+f_an(2).*youtM).*youtM./(4*pi.*modx2).*(1-R^2./modx2);
    
%umod2_an=uout_an.^2+vout_an.^2;
psi=stream2(xoutM,youtM,uoutM,voutM,-2,1);
xy=psi{1};

psi1=stream2(xoutM,youtM,uoutM,voutM,-1,1);
xy1=psi1{1};

psi2=stream2(xoutM,youtM,uoutM,voutM,-0.5,0.5);
xy2=psi2{1};

psi3=stream2(xoutM,youtM,uoutM,voutM,-0.25,0.25);
xy3=psi3{1};

psi4=stream2(xoutM,youtM,uoutM,voutM,-1.0,-1.0);
xy4=psi4{1};

psi5=stream2(xoutM,youtM,uoutM,voutM,-2.0,-2.0);
xy5=psi5{1};

figure
plot(xb, yb, 'k.')
hold on
plot(xy(:,1),xy(:,2),'.')
plot(xy1(:,1),xy1(:,2),'r.')
plot(xy2(:,1),xy2(:,2),'c.')
plot(xy3(:,1),xy3(:,2),'m.')
plot(xy4(:,1),xy4(:,2),'g.')
plot(xy5(:,1),xy5(:,2),'k.')
xlabel('X')
ylabel('Y')
zlabel('Z')
daspect([1 1 1])
grid on
quiver(xoutM,youtM,uoutM,voutM)
contour(xout,yout,((uoutM).^2+voutM.^2),[-1:0.1:6])
axis xy
colorbar
set(gcf,'PaperUnits','centimeters','PaperSize',[20 16],'PaperPosition',[0 0 16 16]);
print(gcf,'-dpdf',['vel_streamlines.pdf'])
 


% figure
% contour(xout,yout,sqrt(uoutM2.^2+voutM2.^2),[-1:0.1:6])
% axis xy
% colorbar

% figure
% imagesc(xout,yout,uoutM-uout_an)
% axis xy
% colorbar
end